knitr::opts_chunk$set(echo = TRUE)
rm(list = ls())
#cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")
#cmdstanr::set_cmdstan_path(path = "C:/Users/pascku/.cmdstan/cmdstan-2.36.0")
library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)
source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE
do_priorsense = FALSE
get_bayesfactor = TRUE
check_models = TRUE #
if (get_bayesfactor) {
stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'BF', 'Rhat', 'ESS')
} else {
stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS')
}
options(
dplyr.print_max = 100,
brms.backend = 'cmdstan',
brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
brms.file_refit = 'on_change',
#brms.file_refit = 'always',
error = function() {
beepr::beep(sound = 5)
if (shutdown) {
system("shutdown /s /t 180")
quit(save = "no", status = 1)
}
}
, es.use_symbols = TRUE
)
####################### Model parameters #######################
iterations = 2000 # 12'000 per chain to achieve 40'000
warmup = 1000 # 2000
# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'
################################################################
suffix = paste0('_as_preregistered_new_', as.character(iterations))df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df
df_double <- prepare_data(df, recode_pushing = TRUE, use_mi = use_mi)[[1]]Constructing scales Re-coding pusing reshaping data (4field) centering data within and between
Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s 0.0000 0.0000 0.0000 0.1649 0.0000 5.0000 275
formula <- bf(
pa_sub ~ 0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
is_A:persuasion_self_cb + is_B:persuasion_self_cb +
is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw |dd| coupleID),
hu = ~ 0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
is_A:persuasion_self_cb + is_B:persuasion_self_cb +
is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_A', dpar = 'hu')
, brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_B', dpar = 'hu')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = hurdle_lognormal()
#)
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_sub_persuasion <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::hurdle_lognormal(),
#family = brms::hurdle_negbinomial(),
#family = brms::hurdle_poisson(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 42,
file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_persuasion", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
check_brms(pa_sub_persuasion, log_pp_check = TRUE)
DHARMa.check_brms.all(pa_sub_persuasion, integer = TRUE, outliers_type = 'bootstrap')
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 2.24 [2.15, 2.34] 1.50 0.45
## is_B 1.06 [1.04, 1.10] 1.03 0.94
## is_A:persuasion_self_cw 1.17 [1.14, 1.21] 1.08 0.86
## is_B:persuasion_self_cw 1.15 [1.12, 1.19] 1.07 0.87
## is_A:persuasion_partner_cw 1.05 [1.03, 1.10] 1.03 0.95
## is_B:persuasion_partner_cw 2.78 [2.66, 2.90] 1.67 0.36
## is_A:persuasion_self_cb 2.83 [2.70, 2.96] 1.68 0.35
## Tolerance 95% CI
## [0.43, 0.46]
## [0.91, 0.96]
## [0.83, 0.88]
## [0.84, 0.90]
## [0.91, 0.97]
## [0.34, 0.38]
## [0.34, 0.37]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 94%
## lognormal 6%
##
## Predicted Distribution of Response
##
## Distribution Probability
## neg. binomial (zero-infl.) 78%
## beta-binomial 16%
## lognormal 3%
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 6, observations = 3742, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.0005344735 0.0030799038
## sample estimates:
## outlier frequency (expected: 0.00170229823623731 )
## 0.001603421
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub_persuasion$data$pa_sub[pa_sub_persuasion$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_sub <- summarize_brms(
pa_sub_persuasion,
stats_to_report = stats_to_report,
rope_range = rope_range_continuous,
hu_rope_range = c(-0.18, 0.18),
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.
## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
## Warning: Bayes factors might not be precise.
## For precise Bayes factors, sampling at least 40,000 posterior
## samples is recommended.
## Warning in summarize_brms(pa_sub_persuasion, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.
| exp(Est.)_hu | SE_hu | 95% CI_hu | pd_hu | ROPE_hu | inside ROPE_hu | BF_hu | BF_Evidence_hu | Rhat_hu | Bulk_ESS_hu | Tail_ESS_hu | exp(Est.)_nonzero | SE_nonzero | 95% CI_nonzero | pd_nonzero | ROPE_nonzero | inside ROPE_nonzero | BF_nonzero | BF_Evidence_nonzero | Rhat_nonzero | Bulk_ESS_nonzero | Tail_ESS_nonzero | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 0.74 | 0.14 | [ 0.51, 1.08] | 0.939 | [0.84, 1.20] | 0.270 | 0.270 | Moderate Evidence for Null | 1.002 | 2269 | 2450 | 46.64*** | 3.64 | [39.95, 53.88] | 1.000 | [0.92, 1.08] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 1427 | 2428 |
| is_A:day | 0.97 | 0.18 | [ 0.67, 1.41] | 0.575 | [0.84, 1.20] | 0.650 | 0.098 | Strong Evidence for Null | 1.001 | 5389 | 2918 | 1.09 | 0.09 | [ 0.92, 1.29] | 0.851 | [0.92, 1.08] | 0.434 | 0.011 | Very Strong Evidence for Null | 1.002 | 4676 | 2918 |
| is_A:persuasion_partner_cb | 0.76 | 0.30 | [ 0.32, 1.62] | 0.765 | [0.84, 1.20] | 0.286 | 0.261 | Moderate Evidence for Null | 1.001 | 1551 | 2445 | 1.03 | 0.17 | [ 0.75, 1.42] | 0.568 | [0.92, 1.08] | 0.365 | 0.017 | Very Strong Evidence for Null | 1.002 | 1518 | 2131 |
| is_A:persuasion_partner_cw | 1.58*** | 0.16 | [ 1.29, 1.97] | 1.000 | [0.84, 1.20] | 0.004 | >100 | Overwhelming Evidence | 1.000 | 2783 | 3038 | 1.02 | 0.03 | [ 0.96, 1.08] | 0.716 | [0.92, 1.08] | 0.981 | 0.008 | Very Strong Evidence for Null | 1.001 | 2593 | 2466 |
| is_A:persuasion_self_cb | 0.81 | 0.38 | [ 0.33, 1.99] | 0.677 | [0.84, 1.20] | 0.274 | 0.268 | Moderate Evidence for Null | 1.000 | 1243 | 2299 | 1.21 | 0.21 | [ 0.84, 1.73] | 0.866 | [0.92, 1.08] | 0.190 | 0.032 | Strong Evidence for Null | 1.002 | 1506 | 2250 |
| is_A:persuasion_self_cw | 1.67*** | 0.16 | [ 1.41, 2.09] | 1.000 | [0.84, 1.20] | 0.000 | >100 | Overwhelming Evidence | 1.002 | 2279 | 1744 | 1.01 | 0.04 | [ 0.94, 1.09] | 0.650 | [0.92, 1.08] | 0.935 | 0.011 | Very Strong Evidence for Null | 1.005 | 1848 | 2443 |
| is_B | 0.92 | 0.18 | [ 0.63, 1.37] | 0.657 | [0.84, 1.20] | 0.589 | 0.083 | Strong Evidence for Null | 1.002 | 2113 | 2747 | 49.88*** | 3.87 | [42.42, 57.98] | 1.000 | [0.92, 1.08] | 0.000 | >100 | Overwhelming Evidence | 1.004 | 1498 | 2350 |
| is_B:day | 0.87 | 0.16 | [ 0.61, 1.27] | 0.776 | [0.84, 1.20] | 0.542 | 0.122 | Moderate Evidence for Null | 1.002 | 5012 | 2585 | 0.88 | 0.07 | [ 0.75, 1.03] | 0.946 | [0.92, 1.08] | 0.246 | 0.027 | Very Strong Evidence for Null | 1.000 | 5060 | 2971 |
| is_B:persuasion_partner_cb | 1.01 | 0.50 | [ 0.38, 2.59] | 0.507 | [0.84, 1.20] | 0.276 | 0.246 | Moderate Evidence for Null | 1.001 | 1212 | 1923 | 1.03 | 0.20 | [ 0.71, 1.48] | 0.562 | [0.92, 1.08] | 0.319 | 0.016 | Very Strong Evidence for Null | 1.003 | 1293 | 1942 |
| is_B:persuasion_partner_cw | 1.48*** | 0.14 | [ 1.24, 1.85] | 1.000 | [0.84, 1.20] | 0.011 | 37.706 | Very Strong Evidence | 1.001 | 3118 | 1930 | 1.03 | 0.03 | [ 0.97, 1.10] | 0.850 | [0.92, 1.08] | 0.929 | 0.015 | Very Strong Evidence for Null | 1.001 | 2987 | 2511 |
| is_B:persuasion_self_cb | 0.70 | 0.30 | [ 0.30, 1.58] | 0.797 | [0.84, 1.20] | 0.234 | 0.282 | Moderate Evidence for Null | 1.001 | 1614 | 2463 | 1.13 | 0.18 | [ 0.81, 1.56] | 0.769 | [0.92, 1.08] | 0.271 | 0.023 | Very Strong Evidence for Null | 1.002 | 1431 | 2204 |
| is_B:persuasion_self_cw | 1.68*** | 0.17 | [ 1.39, 2.12] | 1.000 | [0.84, 1.20] | 0.001 | >100 | Overwhelming Evidence | 1.000 | 2613 | 2076 | 1.03 | 0.03 | [ 0.98, 1.10] | 0.876 | [0.92, 1.08] | 0.932 | 0.015 | Very Strong Evidence for Null | 1.001 | 2680 | 3144 |
| sd(is_A) | 0.94 | 0.13 | [0.73, 1.23] | NA | NA | NA | NA | NA | 1.002 | 1960 | 2814 | 0.33 | 0.05 | [0.25, 0.45] | NA | NA | NA | NA | NA | 1.001 | 1331 | 2454 |
| sd(is_A:persuasion_partner_cw) | 0.42 | 0.11 | [0.21, 0.70] | NA | NA | NA | NA | NA | 1.000 | 1999 | 2257 | 0.06 | 0.03 | [0.01, 0.13] | NA | NA | NA | NA | NA | 1.006 | 1037 | 757 |
| sd(is_A:persuasion_self_cw) | 0.29 | 0.13 | [0.05, 0.59] | NA | NA | NA | NA | NA | 1.003 | 1124 | 1201 | 0.17 | 0.03 | [0.11, 0.23] | NA | NA | NA | NA | NA | 1.004 | 2044 | 2529 |
| sd(is_B) | 1.01 | 0.13 | [0.79, 1.32] | NA | NA | NA | NA | NA | 1.001 | 2288 | 2551 | 0.39 | 0.05 | [0.30, 0.50] | NA | NA | NA | NA | NA | 1.001 | 1665 | 2636 |
| sd(is_B:persuasion_partner_cw) | 0.33 | 0.14 | [0.05, 0.63] | NA | NA | NA | NA | NA | 1.003 | 1118 | 958 | 0.10 | 0.03 | [0.05, 0.17] | NA | NA | NA | NA | NA | 1.003 | 2312 | 2581 |
| sd(is_B:persuasion_self_cw) | 0.38 | 0.11 | [0.19, 0.62] | NA | NA | NA | NA | NA | 1.000 | 2229 | 2918 | 0.09 | 0.03 | [0.05, 0.16] | NA | NA | NA | NA | NA | 1.001 | 2289 | 2531 |
| sigma | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 0.67 | 0.01 | [0.65, 0.70] | NA | NA | NA | NA | NA | 1.001 | 5885 | 2983 |
Hurdle part of the model on the left, non-zero part towards the right side of the table
formula <- bf(
pa_sub ~ 0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw +
is_A:pressure_self_cb + is_B:pressure_self_cb +
is_A:pressure_partner_cb + is_B:pressure_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw |dd| coupleID),
hu = ~ 0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw +
is_A:pressure_self_cb + is_B:pressure_self_cb +
is_A:pressure_partner_cb + is_B:pressure_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_A', dpar = 'hu')
, brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_B', dpar = 'hu')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = hurdle_lognormal()
#)
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_sub_pressure <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::hurdle_lognormal(),
#family = brms::hurdle_negbinomial(),
#family = brms::hurdle_poisson(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 42,
file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_pressure", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
check_brms(pa_sub_pressure, log_pp_check = TRUE)
DHARMa.check_brms.all(pa_sub_pressure, integer = TRUE, outliers_type = 'bootstrap')
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 2.33 [2.23, 2.43] 1.53 0.43
## is_B 1.06 [1.03, 1.10] 1.03 0.95
## is_A:pressure_self_cw 1.02 [1.00, 1.09] 1.01 0.98
## is_B:pressure_self_cw 1.02 [1.00, 1.09] 1.01 0.98
## is_A:pressure_partner_cw 1.02 [1.00, 1.09] 1.01 0.98
## is_B:pressure_partner_cw 3.70 [3.53, 3.88] 1.92 0.27
## is_A:pressure_self_cb 4.13 [3.94, 4.34] 2.03 0.24
## Tolerance 95% CI
## [0.41, 0.45]
## [0.91, 0.97]
## [0.92, 1.00]
## [0.92, 1.00]
## [0.92, 1.00]
## [0.26, 0.28]
## [0.23, 0.25]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 94%
## lognormal 6%
##
## Predicted Distribution of Response
##
## Distribution Probability
## neg. binomial (zero-infl.) 78%
## beta-binomial 16%
## lognormal 3%
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 13 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 6, observations = 3736, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.0005353319 0.0029443255
## sample estimates:
## outlier frequency (expected: 0.00158458244111349 )
## 0.001605996
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub_pressure$data$pa_sub[pa_sub_pressure$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_sub <- summarize_brms(
pa_sub_pressure,
stats_to_report = stats_to_report,
rope_range = rope_range_continuous,
hu_rope_range = c(-0.18, 0.18),
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.
## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
## Warning: Bayes factors might not be precise.
## For precise Bayes factors, sampling at least 40,000 posterior
## samples is recommended.
## Warning in summarize_brms(pa_sub_pressure, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.
| exp(Est.)_hu | SE_hu | 95% CI_hu | pd_hu | ROPE_hu | inside ROPE_hu | BF_hu | BF_Evidence_hu | Rhat_hu | Bulk_ESS_hu | Tail_ESS_hu | exp(Est.)_nonzero | SE_nonzero | 95% CI_nonzero | pd_nonzero | ROPE_nonzero | inside ROPE_nonzero | BF_nonzero | BF_Evidence_nonzero | Rhat_nonzero | Bulk_ESS_nonzero | Tail_ESS_nonzero | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 0.86 | 0.16 | [ 0.62, 1.22] | 0.778 | [0.84, 1.20] | 0.543 | 0.093 | Strong Evidence for Null | 1.002 | 1619 | 2259 | 48.30*** | 3.53 | [41.60, 55.61] | 1.000 | [0.92, 1.08] | 0.000 | >100 | Overwhelming Evidence | 1.001 | 1078 | 1836 |
| is_A:day | 0.68* | 0.11 | [ 0.49, 0.94] | 0.992 | [0.84, 1.20] | 0.102 | 1.264 | Weak Evidence | 1.000 | 3627 | 2877 | 1.09 | 0.09 | [ 0.93, 1.27] | 0.835 | [0.92, 1.08] | 0.455 | 0.011 | Very Strong Evidence for Null | 1.000 | 5148 | 3217 |
| is_A:pressure_partner_cb | 0.63 | 0.38 | [ 0.19, 2.19] | 0.768 | [0.84, 1.20] | 0.170 | 0.441 | Weak Evidence for Null | 1.002 | 1387 | 2231 | 0.94 | 0.27 | [ 0.54, 1.65] | 0.581 | [0.92, 1.08] | 0.210 | 0.017 | Very Strong Evidence for Null | 1.004 | 1193 | 2193 |
| is_A:pressure_partner_cw | 2.43*** | 0.57 | [ 1.58, 4.20] | 1.000 | [0.84, 1.20] | 0.000 | >100 | Overwhelming Evidence | 1.001 | 3014 | 2270 | 0.94 | 0.05 | [ 0.84, 1.05] | 0.882 | [0.92, 1.08] | 0.574 | 0.013 | Very Strong Evidence for Null | 1.000 | 3919 | 2837 |
| is_A:pressure_self_cb | 0.56 | 0.43 | [ 0.13, 2.53] | 0.774 | [0.84, 1.20] | 0.143 | 0.497 | Weak Evidence for Null | 1.001 | 1507 | 2224 | 1.56 | 0.52 | [ 0.79, 3.18] | 0.907 | [0.92, 1.08] | 0.071 | 0.040 | Strong Evidence for Null | 1.005 | 1500 | 2130 |
| is_A:pressure_self_cw | 1.82* | 0.59 | [ 1.03, 4.33] | 0.978 | [0.84, 1.20] | 0.061 | 1.476 | Weak Evidence | 1.001 | 2051 | 1761 | 0.85* | 0.06 | [ 0.72, 0.99] | 0.980 | [0.92, 1.08] | 0.135 | 0.079 | Strong Evidence for Null | 1.001 | 2715 | 2077 |
| is_B | 1.14 | 0.22 | [ 0.80, 1.65] | 0.769 | [0.84, 1.20] | 0.551 | 0.095 | Strong Evidence for Null | 1.002 | 1294 | 1946 | 52.19*** | 4.31 | [44.52, 61.49] | 1.000 | [0.92, 1.08] | 0.000 | >100 | Overwhelming Evidence | 1.003 | 875 | 1955 |
| is_B:day | 0.56*** | 0.10 | [ 0.39, 0.79] | 1.000 | [0.84, 1.20] | 0.011 | 17.113 | Strong Evidence | 1.000 | 4000 | 3081 | 0.83* | 0.07 | [ 0.72, 0.97] | 0.990 | [0.92, 1.08] | 0.090 | 0.109 | Moderate Evidence for Null | 1.000 | 4854 | 3016 |
| is_B:pressure_partner_cb | 1.15 | 0.87 | [ 0.26, 5.29] | 0.568 | [0.84, 1.20] | 0.186 | 0.357 | Weak Evidence for Null | 1.001 | 1340 | 1906 | 1.09 | 0.44 | [ 0.48, 2.50] | 0.586 | [0.92, 1.08] | 0.156 | 0.020 | Very Strong Evidence for Null | 1.003 | 1237 | 1738 |
| is_B:pressure_partner_cw | 3.33* | 2.03 | [ 1.23, 16.18] | 0.989 | [0.84, 1.20] | 0.017 | 4.570 | Moderate Evidence | 1.001 | 1929 | 2039 | 0.99 | 0.07 | [ 0.82, 1.13] | 0.567 | [0.92, 1.08] | 0.718 | 0.007 | Very Strong Evidence for Null | 1.000 | 2374 | 2131 |
| is_B:pressure_self_cb | 0.48 | 0.29 | [ 0.13, 1.54] | 0.893 | [0.84, 1.20] | 0.116 | 0.615 | Weak Evidence for Null | 1.001 | 1311 | 2034 | 1.08 | 0.34 | [ 0.55, 2.03] | 0.591 | [0.92, 1.08] | 0.189 | 0.019 | Very Strong Evidence for Null | 1.001 | 1278 | 2026 |
| is_B:pressure_self_cw | 1.29 | 0.35 | [ 0.65, 2.65] | 0.828 | [0.84, 1.20] | 0.316 | 0.229 | Moderate Evidence for Null | 1.001 | 2208 | 1848 | 0.94 | 0.06 | [ 0.82, 1.08] | 0.834 | [0.92, 1.08] | 0.548 | 0.012 | Very Strong Evidence for Null | 1.002 | 2727 | 2429 |
| sd(is_A) | 0.88 | 0.12 | [0.69, 1.15] | NA | NA | NA | NA | NA | 1.001 | 1982 | 2636 | 0.32 | 0.05 | [0.24, 0.43] | NA | NA | NA | NA | NA | 1.004 | 1243 | 2039 |
| sd(is_A:pressure_partner_cw) | 0.30 | 0.28 | [0.01, 1.22] | NA | NA | NA | NA | NA | 1.001 | 1672 | 1800 | 0.06 | 0.05 | [0.00, 0.21] | NA | NA | NA | NA | NA | 1.000 | 2120 | 2016 |
| sd(is_A:pressure_self_cw) | 0.70 | 0.51 | [0.06, 2.24] | NA | NA | NA | NA | NA | 1.004 | 1258 | 1332 | 0.10 | 0.08 | [0.00, 0.37] | NA | NA | NA | NA | NA | 1.001 | 1512 | 1566 |
| sd(is_B) | 0.91 | 0.12 | [0.70, 1.19] | NA | NA | NA | NA | NA | 1.002 | 1828 | 2298 | 0.40 | 0.06 | [0.31, 0.52] | NA | NA | NA | NA | NA | 1.005 | 1150 | 2250 |
| sd(is_B:pressure_partner_cw) | 1.67 | 0.87 | [0.36, 3.91] | NA | NA | NA | NA | NA | 1.002 | 1273 | 960 | 0.09 | 0.08 | [0.00, 0.37] | NA | NA | NA | NA | NA | 1.001 | 1324 | 1668 |
| sd(is_B:pressure_self_cw) | 0.75 | 0.49 | [0.06, 2.09] | NA | NA | NA | NA | NA | 1.005 | 1011 | 980 | 0.09 | 0.07 | [0.01, 0.28] | NA | NA | NA | NA | NA | 1.001 | 1633 | 1710 |
| sigma | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 0.69 | 0.01 | [0.67, 0.72] | NA | NA | NA | NA | NA | 1.000 | 4519 | 3405 |
Hurdle part of the model on the left, non-zero part towards the right side of the table
formula <- bf(
pa_sub ~ 0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw +
is_A:pushing_self_cb + is_B:pushing_self_cb +
is_A:pushing_partner_cb + is_B:pushing_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
barriers_self_cw + barriers_self_cb +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID),
hu = ~ 0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw +
is_A:pushing_self_cb + is_B:pushing_self_cb +
is_A:pushing_partner_cb + is_B:pushing_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_A', dpar = 'hu')
, brms::set_prior("normal(0.5, 2.5)", class = "b", coef = 'is_B', dpar = 'hu')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = hurdle_lognormal()
#)
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_sub_pushing <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::hurdle_lognormal(),
#family = brms::hurdle_negbinomial(),
#family = brms::hurdle_poisson(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 42,
file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_pushing", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
check_brms(pa_sub_pushing, log_pp_check = TRUE)
DHARMa.check_brms.all(pa_sub_pushing, integer = TRUE, outliers_type = 'bootstrap')
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.84 [1.76, 1.93] 1.36 0.54
## is_B 1.09 [1.06, 1.14] 1.04 0.92
## barriers_self_cw 1.05 [1.02, 1.10] 1.02 0.95
## barriers_self_cb 1.04 [1.02, 1.10] 1.02 0.96
## is_A:pushing_self_cw 1.08 [1.05, 1.13] 1.04 0.92
## is_B:pushing_self_cw 1.23 [1.19, 1.28] 1.11 0.81
## is_A:pushing_partner_cw 1.31 [1.26, 1.37] 1.15 0.76
## is_B:pushing_partner_cw 1.24 [1.20, 1.29] 1.11 0.81
## is_A:pushing_self_cb 1.28 [1.23, 1.33] 1.13 0.78
## Tolerance 95% CI
## [0.52, 0.57]
## [0.88, 0.94]
## [0.91, 0.98]
## [0.91, 0.98]
## [0.88, 0.95]
## [0.78, 0.84]
## [0.73, 0.79]
## [0.77, 0.84]
## [0.75, 0.81]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 91%
## lognormal 9%
##
## Predicted Distribution of Response
##
## Distribution Probability
## neg. binomial (zero-infl.) 84%
## beta-binomial 9%
## lognormal 6%
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 10 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 3, observations = 1776, p-value =
## 0.94
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.003673986
## sample estimates:
## outlier frequency (expected: 0.00148085585585586 )
## 0.001689189
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub_pushing$data$pa_sub[pa_sub_pushing$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_sub <- summarize_brms(
pa_sub_pushing,
stats_to_report = stats_to_report,
rope_range = rope_range_continuous,
hu_rope_range = c(-0.18, 0.18),
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.
## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
## Warning: Bayes factors might not be precise.
## For precise Bayes factors, sampling at least 40,000 posterior
## samples is recommended.
## Warning in summarize_brms(pa_sub_pushing, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.
| exp(Est.)_hu | SE_hu | 95% CI_hu | pd_hu | ROPE_hu | inside ROPE_hu | BF_hu | BF_Evidence_hu | Rhat_hu | Bulk_ESS_hu | Tail_ESS_hu | exp(Est.)_nonzero | SE_nonzero | 95% CI_nonzero | pd_nonzero | ROPE_nonzero | inside ROPE_nonzero | BF_nonzero | BF_Evidence_nonzero | Rhat_nonzero | Bulk_ESS_nonzero | Tail_ESS_nonzero | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| barriers_self_cb | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 0.97 | 0.13 | [ 0.75, 1.27] | 0.605 | [0.93, 1.08] | 0.416 | 0.019 | Very Strong Evidence for Null | 1.001 | 1764 | 2409 |
| barriers_self_cw | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 1.07* | 0.03 | [ 1.01, 1.13] | 0.985 | [0.93, 1.08] | 0.651 | 0.093 | Strong Evidence for Null | 1.001 | 7429 | 3025 |
| is_A | 2.61*** | 0.68 | [ 1.57, 4.45] | 1.000 | [0.84, 1.20] | 0.002 | 49.120 | Very Strong Evidence | 1.001 | 2799 | 3158 | 53.29*** | 4.62 | [44.75, 63.15] | 1.000 | [0.93, 1.08] | 0.000 | >100 | Overwhelming Evidence | 1.001 | 1898 | 2276 |
| is_A:day | 0.82 | 0.25 | [ 0.46, 1.45] | 0.748 | [0.84, 1.20] | 0.375 | 0.182 | Moderate Evidence for Null | 0.999 | 5470 | 3172 | 1.08 | 0.11 | [ 0.89, 1.31] | 0.795 | [0.93, 1.08] | 0.433 | 0.011 | Very Strong Evidence for Null | 1.000 | 5968 | 2947 |
| is_A:pushing_partner_cb | 0.91 | 0.89 | [ 0.11, 6.51] | 0.537 | [0.84, 1.20] | 0.146 | 0.463 | Weak Evidence for Null | 1.001 | 2860 | 3041 | 1.14 | 0.41 | [ 0.55, 2.32] | 0.638 | [0.93, 1.08] | 0.149 | 0.021 | Very Strong Evidence for Null | 1.002 | 1511 | 2291 |
| is_A:pushing_partner_cw | 1.96*** | 0.44 | [ 1.30, 3.48] | 1.000 | [0.84, 1.20] | 0.010 | 19.295 | Strong Evidence | 1.000 | 3298 | 2482 | 0.96 | 0.03 | [ 0.89, 1.03] | 0.890 | [0.93, 1.08] | 0.833 | 0.019 | Very Strong Evidence for Null | 1.002 | 4654 | 2884 |
| is_A:pushing_self_cb | 0.39 | 0.34 | [ 0.07, 2.33] | 0.861 | [0.84, 1.20] | 0.088 | 0.808 | Weak Evidence for Null | 1.000 | 3309 | 2824 | 1.28 | 0.41 | [ 0.69, 2.46] | 0.777 | [0.93, 1.08] | 0.145 | 0.025 | Very Strong Evidence for Null | 1.001 | 1828 | 2493 |
| is_A:pushing_self_cw | 1.93 | 0.79 | [ 0.94, 5.66] | 0.964 | [0.84, 1.20] | 0.074 | 1.112 | Weak Evidence | 1.001 | 2503 | 2226 | 1.00 | 0.06 | [ 0.87, 1.12] | 0.530 | [0.93, 1.08] | 0.798 | 0.012 | Very Strong Evidence for Null | 1.001 | 2157 | 2444 |
| is_B | 4.00*** | 1.26 | [ 2.24, 7.32] | 1.000 | [0.84, 1.20] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 3081 | 2528 | 52.70*** | 4.59 | [44.27, 62.31] | 1.000 | [0.93, 1.08] | 0.000 | >100 | Overwhelming Evidence | 1.001 | 1867 | 2678 |
| is_B:day | 0.64 | 0.20 | [ 0.36, 1.16] | 0.920 | [0.84, 1.20] | 0.186 | 0.390 | Weak Evidence for Null | 1.002 | 6963 | 2937 | 0.94 | 0.09 | [ 0.79, 1.13] | 0.759 | [0.93, 1.08] | 0.492 | 0.010 | Very Strong Evidence for Null | 1.001 | 8205 | 2886 |
| is_B:pushing_partner_cb | 0.79 | 0.78 | [ 0.10, 5.57] | 0.592 | [0.84, 1.20] | 0.142 | 0.472 | Weak Evidence for Null | 1.001 | 2790 | 2898 | 1.03 | 0.35 | [ 0.52, 2.05] | 0.533 | [0.93, 1.08] | 0.177 | 0.022 | Very Strong Evidence for Null | 1.002 | 2021 | 2507 |
| is_B:pushing_partner_cw | 2.32* | 1.13 | [ 1.04, 8.50] | 0.979 | [0.84, 1.20] | 0.040 | 1.872 | Weak Evidence | 1.001 | 2151 | 2430 | 1.02 | 0.06 | [ 0.91, 1.14] | 0.660 | [0.93, 1.08] | 0.778 | 0.013 | Very Strong Evidence for Null | 1.000 | 3161 | 2928 |
| is_B:pushing_self_cb | 0.40 | 0.41 | [ 0.05, 3.16] | 0.807 | [0.84, 1.20] | 0.084 | 0.881 | Weak Evidence for Null | 1.001 | 3290 | 3081 | 1.47 | 0.59 | [ 0.67, 3.28] | 0.847 | [0.93, 1.08] | 0.095 | 0.037 | Strong Evidence for Null | 1.000 | 1905 | 2362 |
| is_B:pushing_self_cw | 1.16 | 0.25 | [ 0.72, 1.90] | 0.745 | [0.84, 1.20] | 0.477 | 0.158 | Moderate Evidence for Null | 1.001 | 2795 | 2632 | 0.94 | 0.03 | [ 0.88, 1.01] | 0.943 | [0.93, 1.08] | 0.708 | 0.028 | Very Strong Evidence for Null | 1.000 | 4564 | 3152 |
| sd(is_A) | 1.17 | 0.20 | [0.86, 1.65] | NA | NA | NA | NA | NA | 1.001 | 2876 | 3069 | 0.36 | 0.05 | [0.26, 0.50] | NA | NA | NA | NA | NA | 1.001 | 1765 | 2700 |
| sd(is_A:pushing_partner_cw) | 0.66 | 0.33 | [0.10, 1.48] | NA | NA | NA | NA | NA | 1.000 | 1570 | 1332 | 0.04 | 0.04 | [0.00, 0.14] | NA | NA | NA | NA | NA | 1.003 | 1972 | 2212 |
| sd(is_A:pushing_self_cw) | 1.54 | 0.59 | [0.51, 3.12] | NA | NA | NA | NA | NA | 1.002 | 1619 | 1623 | 0.18 | 0.06 | [0.06, 0.33] | NA | NA | NA | NA | NA | 1.001 | 1388 | 1087 |
| sd(is_B) | 1.38 | 0.23 | [0.99, 1.92] | NA | NA | NA | NA | NA | 1.000 | 2820 | 2960 | 0.39 | 0.06 | [0.29, 0.54] | NA | NA | NA | NA | NA | 1.001 | 1811 | 2685 |
| sd(is_B:pushing_partner_cw) | 1.74 | 0.84 | [0.24, 3.69] | NA | NA | NA | NA | NA | 1.002 | 1076 | 922 | 0.17 | 0.07 | [0.03, 0.31] | NA | NA | NA | NA | NA | 1.001 | 1494 | 1330 |
| sd(is_B:pushing_self_cw) | 0.85 | 0.32 | [0.32, 1.74] | NA | NA | NA | NA | NA | 1.001 | 1714 | 1917 | 0.05 | 0.04 | [0.00, 0.15] | NA | NA | NA | NA | NA | 1.001 | 1790 | 2025 |
| sigma | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 0.67 | 0.01 | [0.64, 0.70] | NA | NA | NA | NA | NA | 1.001 | 4994 | 2914 |
Hurdle part of the model on the left, non-zero part towards the right side of the table
## [1] 5.75 971.25
formula <- bf(
pa_obj ~ 0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
is_A:persuasion_self_cb + is_B:persuasion_self_cb +
is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
is_A:weartime_self_cw + is_B:weartime_self_cw +
is_A:weartime_self_cb + is_B:weartime_self_cb +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = lognormal()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_obj_log_persuasion <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = lognormal(),
#control = list(adapt_delta = 0.95),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_persuasion", suffix))
)## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
check_brms(pa_obj_log_persuasion, log_pp_check = TRUE)
DHARMa.check_brms.all(pa_obj_log_persuasion, integer = TRUE, outliers_type = 'bootstrap')
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.94 [1.86, 2.04] 1.39 0.51
## is_B 1.09 [1.06, 1.14] 1.04 0.92
## is_A:persuasion_self_cw 1.17 [1.13, 1.21] 1.08 0.86
## is_B:persuasion_self_cw 1.12 [1.08, 1.16] 1.06 0.90
## is_A:persuasion_partner_cw 1.05 [1.02, 1.10] 1.02 0.95
## is_B:persuasion_partner_cw 3.62 [3.43, 3.82] 1.90 0.28
## is_A:persuasion_self_cb 2.84 [2.70, 2.99] 1.68 0.35
## is_B:persuasion_self_cb 3.82 [3.62, 4.04] 1.95 0.26
## is_A:persuasion_partner_cb 2.87 [2.72, 3.02] 1.69 0.35
## Tolerance 95% CI
## [0.49, 0.54]
## [0.88, 0.94]
## [0.82, 0.88]
## [0.86, 0.92]
## [0.91, 0.98]
## [0.26, 0.29]
## [0.33, 0.37]
## [0.25, 0.28]
## [0.33, 0.37]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 88%
## lognormal 9%
## weibull 3%
##
## Predicted Distribution of Response
##
## Distribution Probability
## lognormal 72%
## tweedie 16%
## neg. binomial (zero-infl.) 9%
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 31, observations = 3343, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.001495663 0.004786120
## sample estimates:
## outlier frequency (expected: 0.00320370924319474 )
## 0.009273108
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log_persuasion$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_obj <- summarize_brms(
pa_obj_log_persuasion,
stats_to_report = stats_to_report,
rope_range = rope_range_log,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
## Warning: Bayes factors might not be precise.
## For precise Bayes factors, sampling at least 40,000 posterior
## samples is recommended.
## Warning in summarize_brms(pa_obj_log_persuasion, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.
| exp(Est.) | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 101.71*** | 6.22 | [ 89.97, 115.11] | 1.000 | [0.94, 1.07] | 0.000 | >100 | Overwhelming Evidence | 1.003 | 951 | 1434 |
| is_B | 136.08*** | 7.27 | [122.13, 151.56] | 1.000 | [0.94, 1.07] | 0.000 | >100 | Overwhelming Evidence | 1.005 | 1389 | 1918 |
| is_A:persuasion_self_cw | 1.05* | 0.02 | [ 1.01, 1.09] | 0.987 | [0.94, 1.07] | 0.843 | 0.076 | Strong Evidence for Null | 1.000 | 3928 | 3337 |
| is_B:persuasion_self_cw | 1.01 | 0.02 | [ 0.97, 1.05] | 0.655 | [0.94, 1.07] | 0.997 | 0.006 | Very Strong Evidence for Null | 1.001 | 3630 | 3130 |
| is_A:persuasion_partner_cw | 1.02 | 0.02 | [ 0.98, 1.06] | 0.796 | [0.94, 1.07] | 0.987 | 0.008 | Very Strong Evidence for Null | 0.999 | 3459 | 3176 |
| is_B:persuasion_partner_cw | 1.03 | 0.02 | [ 0.99, 1.07] | 0.913 | [0.94, 1.07] | 0.964 | 0.014 | Very Strong Evidence for Null | 1.000 | 3753 | 3537 |
| is_A:persuasion_self_cb | 1.25 | 0.22 | [ 0.87, 1.77] | 0.899 | [0.94, 1.07] | 0.128 | 0.036 | Strong Evidence for Null | 1.004 | 994 | 1743 |
| is_B:persuasion_self_cb | 0.89 | 0.11 | [ 0.69, 1.13] | 0.837 | [0.94, 1.07] | 0.264 | 0.018 | Very Strong Evidence for Null | 1.003 | 1266 | 1932 |
| is_A:persuasion_partner_cb | 0.89 | 0.14 | [ 0.67, 1.21] | 0.785 | [0.94, 1.07] | 0.244 | 0.021 | Very Strong Evidence for Null | 1.004 | 1106 | 2003 |
| is_B:persuasion_partner_cb | 1.29 | 0.19 | [ 0.97, 1.76] | 0.960 | [0.94, 1.07] | 0.086 | 0.059 | Strong Evidence for Null | 1.003 | 1153 | 2022 |
| is_A:day | 0.98 | 0.05 | [ 0.90, 1.08] | 0.626 | [0.94, 1.07] | 0.805 | 0.004 | Very Strong Evidence for Null | 1.001 | 9420 | 2738 |
| is_B:day | 0.95 | 0.05 | [ 0.86, 1.04] | 0.856 | [0.94, 1.07] | 0.608 | 0.006 | Very Strong Evidence for Null | 1.003 | 8586 | 2175 |
| is_A:weartime_self_cw | 1.00*** | 0.00 | [ 1.00, 1.00] | 1.000 | [0.94, 1.07] | 1.000 | 3.072 | Moderate Evidence | 1.000 | 9905 | 3242 |
| is_B:weartime_self_cw | 1.00* | 0.00 | [ 1.00, 1.00] | 0.986 | [0.94, 1.07] | 1.000 | 0.042 | Strong Evidence for Null | 1.002 | 7670 | 2893 |
| is_A:weartime_self_cb | 1.00 | 0.00 | [ 1.00, 1.00] | 0.702 | [0.94, 1.07] | 1.000 | 0.013 | Very Strong Evidence for Null | 1.002 | 1124 | 1965 |
| is_B:weartime_self_cb | 1.00 | 0.00 | [ 1.00, 1.00] | 0.575 | [0.94, 1.07] | 1.000 | 0.011 | Very Strong Evidence for Null | 1.001 | 1304 | 1992 |
| sd(is_A) | 0.34 | 0.04 | [0.27, 0.44] | NA | NA | NA | NA | NA | 1.001 | 1313 | 2072 |
| sd(is_B) | 0.29 | 0.04 | [0.22, 0.38] | NA | NA | NA | NA | NA | 1.001 | 1472 | 2002 |
| sd(is_A:persuasion_self_cw) | 0.06 | 0.02 | [0.01, 0.11] | NA | NA | NA | NA | NA | 1.006 | 851 | 503 |
| sd(is_B:persuasion_self_cw) | 0.06 | 0.02 | [0.02, 0.11] | NA | NA | NA | NA | NA | 1.002 | 1911 | 1806 |
| sd(is_A:persuasion_partner_cw) | 0.07 | 0.02 | [0.02, 0.13] | NA | NA | NA | NA | NA | 1.001 | 1384 | 1295 |
| sd(is_B:persuasion_partner_cw) | 0.07 | 0.03 | [0.01, 0.13] | NA | NA | NA | NA | NA | 1.002 | 1208 | 1198 |
| sigma | 0.55 | 0.01 | [0.53, 0.56] | NA | NA | NA | NA | NA | 1.000 | 8451 | 3220 |
formula <- bf(
pa_obj ~ 0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw +
is_A:pressure_self_cb + is_B:pressure_self_cb +
is_A:pressure_partner_cb + is_B:pressure_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
is_A:weartime_self_cw + is_B:weartime_self_cw +
is_A:weartime_self_cb + is_B:weartime_self_cb +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = lognormal()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_obj_log_pressure <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = lognormal(),
#control = list(adapt_delta = 0.95),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_pressure", suffix))
)## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
check_brms(pa_obj_log_pressure, log_pp_check = TRUE)
DHARMa.check_brms.all(pa_obj_log_pressure, integer = TRUE, outliers_type = 'bootstrap')
}## Model has no intercept. VIFs may not be sensible.
## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model without
## interaction terms.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.90 [ 1.82, 1.99] 1.38 0.53
## is_B 1.03 [ 1.01, 1.09] 1.01 0.97
## is_A:pressure_self_cw 1.02 [ 1.00, 1.11] 1.01 0.98
## is_B:pressure_self_cw 1.03 [ 1.01, 1.09] 1.01 0.97
## is_A:pressure_partner_cw 1.02 [ 1.01, 1.10] 1.01 0.98
## Tolerance 95% CI
## [0.50, 0.55]
## [0.92, 0.99]
## [0.90, 1.00]
## [0.92, 0.99]
## [0.91, 0.99]
##
## Moderate Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A:pressure_self_cb 8.76 [ 8.26, 9.29] 2.96 0.11
## is_A:pressure_partner_cb 9.02 [ 8.51, 9.57] 3.00 0.11
## Tolerance 95% CI
## [0.11, 0.12]
## [0.10, 0.12]
##
## High Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_B:pressure_partner_cw 10.59 [ 9.98, 11.24] 3.25 0.09
## is_B:pressure_self_cb 10.72 [10.11, 11.38] 3.27 0.09
## Tolerance 95% CI
## [0.09, 0.10]
## [0.09, 0.10]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 88%
## lognormal 9%
## weibull 3%
##
## Predicted Distribution of Response
##
## Distribution Probability
## lognormal 72%
## tweedie 16%
## neg. binomial (zero-infl.) 9%
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 3 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 29, observations = 3337, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.001498352 0.005094396
## sample estimates:
## outlier frequency (expected: 0.00310758166017381 )
## 0.008690441
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log_pressure$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_obj <- summarize_brms(
pa_obj_log_pressure,
stats_to_report = stats_to_report,
rope_range = rope_range_log,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.
## Model has interaction terms. VIFs might be inflated.
## You may check multicollinearity among predictors of a model without
## interaction terms.
## Warning in compute_rope(range = rope_range): Collinearity detected. Some
## VIFs are > 10. This may invalidate ROPE inferences!
Low Correlation
Term VIF VIF 95% CI Increased SE Tolerance
is_A 1.90 [ 1.82, 1.99] 1.38 0.53
is_B 1.03 [ 1.01, 1.09] 1.01 0.97
is_A:pressure_self_cw 1.02 [ 1.00, 1.11] 1.01 0.98
is_B:pressure_self_cw 1.03 [ 1.01, 1.09] 1.01 0.97
is_A:pressure_partner_cw 1.02 [ 1.01, 1.10] 1.01 0.98 Tolerance 95% CI [0.50, 0.55] [0.92, 0.99] [0.90, 1.00] [0.92, 0.99] [0.91, 0.99]
Moderate Correlation
Term VIF VIF 95% CI Increased SE Tolerance
is_A:pressure_self_cb 8.76 [ 8.26, 9.29] 2.96 0.11
is_A:pressure_partner_cb 9.02 [ 8.51, 9.57] 3.00 0.11 Tolerance 95% CI [0.11, 0.12] [0.10, 0.12]
High Correlation
Term VIF VIF 95% CI Increased SE Tolerance
is_B:pressure_partner_cw 10.59 [ 9.98, 11.24] 3.25 0.09 is_B:pressure_self_cb 10.72 [10.11, 11.38] 3.27 0.09 Tolerance 95% CI [0.09, 0.10] [0.09, 0.10]
## Sampling priors, please wait...
## Warning: Bayes factors might not be precise.
## For precise Bayes factors, sampling at least 40,000 posterior
## samples is recommended.
## Warning in summarize_brms(pa_obj_log_pressure, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.
| exp(Est.) | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 101.82*** | 6.29 | [ 90.01, 115.55] | 1.000 | [0.94, 1.07] | 0.000 | >100 | Overwhelming Evidence | 1.001 | 760 | 1414 |
| is_B | 137.46*** | 7.80 | [123.29, 153.89] | 1.000 | [0.94, 1.07] | 0.000 | >100 | Overwhelming Evidence | 1.001 | 939 | 2113 |
| is_A:pressure_self_cw | 0.99 | 0.05 | [ 0.87, 1.09] | 0.607 | [0.94, 1.07] | 0.760 | 0.005 | Very Strong Evidence for Null | 1.001 | 3475 | 2940 |
| is_B:pressure_self_cw | 0.99 | 0.04 | [ 0.91, 1.08] | 0.609 | [0.94, 1.07] | 0.853 | 0.005 | Very Strong Evidence for Null | 0.999 | 5396 | 3072 |
| is_A:pressure_partner_cw | 1.01 | 0.04 | [ 0.93, 1.09] | 0.600 | [0.94, 1.07] | 0.898 | 0.004 | Very Strong Evidence for Null | 1.000 | 5928 | 3035 |
| is_B:pressure_partner_cw | 1.02 | 0.06 | [ 0.92, 1.16] | 0.665 | [0.94, 1.07] | 0.713 | 0.006 | Very Strong Evidence for Null | 1.001 | 3788 | 2687 |
| is_A:pressure_self_cb | 0.85 | 0.29 | [ 0.43, 1.66] | 0.690 | [0.94, 1.07] | 0.132 | 0.017 | Very Strong Evidence for Null | 1.004 | 1534 | 1963 |
| is_B:pressure_self_cb | 0.93 | 0.23 | [ 0.58, 1.53] | 0.618 | [0.94, 1.07] | 0.202 | 0.012 | Very Strong Evidence for Null | 1.003 | 1593 | 2174 |
| is_A:pressure_partner_cb | 1.24 | 0.34 | [ 0.73, 2.15] | 0.791 | [0.94, 1.07] | 0.141 | 0.021 | Very Strong Evidence for Null | 1.005 | 1437 | 2207 |
| is_B:pressure_partner_cb | 1.21 | 0.38 | [ 0.64, 2.23] | 0.726 | [0.94, 1.07] | 0.123 | 0.016 | Very Strong Evidence for Null | 1.004 | 1678 | 2434 |
| is_A:day | 0.95 | 0.04 | [ 0.87, 1.04] | 0.851 | [0.94, 1.07] | 0.642 | 0.006 | Very Strong Evidence for Null | 1.002 | 8015 | 3016 |
| is_B:day | 0.92 | 0.04 | [ 0.84, 1.01] | 0.964 | [0.94, 1.07] | 0.352 | 0.019 | Very Strong Evidence for Null | 1.003 | 9930 | 2577 |
| is_A:weartime_self_cw | 1.00*** | 0.00 | [ 1.00, 1.00] | 1.000 | [0.94, 1.07] | 1.000 | 3.889 | Moderate Evidence | 1.000 | 7074 | 2735 |
| is_B:weartime_self_cw | 1.00 | 0.00 | [ 1.00, 1.00] | 0.971 | [0.94, 1.07] | 1.000 | 0.026 | Very Strong Evidence for Null | 1.001 | 8048 | 2482 |
| is_A:weartime_self_cb | 1.00 | 0.00 | [ 1.00, 1.00] | 0.626 | [0.94, 1.07] | 1.000 | 0.011 | Very Strong Evidence for Null | 1.003 | 914 | 1459 |
| is_B:weartime_self_cb | 1.00 | 0.00 | [ 1.00, 1.00] | 0.718 | [0.94, 1.07] | 1.000 | 0.012 | Very Strong Evidence for Null | 1.009 | 1187 | 2234 |
| sd(is_A) | 0.34 | 0.04 | [0.27, 0.44] | NA | NA | NA | NA | NA | 1.003 | 1060 | 1993 |
| sd(is_B) | 0.29 | 0.04 | [0.23, 0.39] | NA | NA | NA | NA | NA | 1.002 | 1455 | 2240 |
| sd(is_A:pressure_self_cw) | 0.08 | 0.07 | [0.00, 0.25] | NA | NA | NA | NA | NA | 1.000 | 1310 | 1811 |
| sd(is_B:pressure_self_cw) | 0.05 | 0.04 | [0.00, 0.18] | NA | NA | NA | NA | NA | 1.001 | 1935 | 1832 |
| sd(is_A:pressure_partner_cw) | 0.04 | 0.03 | [0.00, 0.13] | NA | NA | NA | NA | NA | 1.002 | 2463 | 1819 |
| sd(is_B:pressure_partner_cw) | 0.07 | 0.06 | [0.00, 0.26] | NA | NA | NA | NA | NA | 1.000 | 1722 | 1765 |
| sigma | 0.56 | 0.01 | [0.54, 0.57] | NA | NA | NA | NA | NA | 1.001 | 7882 | 2946 |
formula <- bf(
pa_obj ~ 0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw +
is_A:pushing_self_cb + is_B:pushing_self_cb +
is_A:pushing_partner_cb + is_B:pushing_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
barriers_self_cw + barriers_self_cb +
is_A:day + is_B:day +
is_A:weartime_self_cw + is_B:weartime_self_cw +
is_A:weartime_self_cb + is_B:weartime_self_cb +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = lognormal()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_obj_log_pushing <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = lognormal(),
#control = list(adapt_delta = 0.95),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_pushing", suffix))
)## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
check_brms(pa_obj_log_pushing, log_pp_check = TRUE)
DHARMa.check_brms.all(pa_obj_log_pushing, integer = TRUE, outliers_type = 'bootstrap')
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 2.03 [1.92, 2.17] 1.43 0.49
## is_B 1.02 [1.00, 1.14] 1.01 0.98
## barriers_self_cw 1.09 [1.05, 1.15] 1.04 0.92
## barriers_self_cb 1.07 [1.03, 1.13] 1.03 0.94
## is_A:pushing_self_cw 1.06 [1.03, 1.13] 1.03 0.94
## is_B:pushing_self_cw 1.05 [1.02, 1.12] 1.02 0.96
## is_A:pushing_partner_cw 1.10 [1.06, 1.16] 1.05 0.91
## is_B:pushing_partner_cw 1.91 [1.80, 2.03] 1.38 0.52
## is_A:pushing_self_cb 1.82 [1.72, 1.94] 1.35 0.55
## is_B:pushing_self_cb 1.93 [1.82, 2.05] 1.39 0.52
## is_A:pushing_partner_cb 1.80 [1.70, 1.91] 1.34 0.56
## Tolerance 95% CI
## [0.46, 0.52]
## [0.88, 1.00]
## [0.87, 0.95]
## [0.88, 0.97]
## [0.89, 0.97]
## [0.89, 0.98]
## [0.86, 0.94]
## [0.49, 0.56]
## [0.52, 0.58]
## [0.49, 0.55]
## [0.52, 0.59]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 88%
## lognormal 9%
## weibull 3%
##
## Predicted Distribution of Response
##
## Distribution Probability
## lognormal 72%
## tweedie 16%
## neg. binomial (zero-infl.) 9%
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 3 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 14, observations = 1594, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.0006273526 0.0062735257
## sample estimates:
## outlier frequency (expected: 0.00330614805520703 )
## 0.008782936
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log_pushing$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_obj <- summarize_brms(
pa_obj_log_pushing,
stats_to_report = stats_to_report,
rope_range = rope_range_log,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
## Warning: Bayes factors might not be precise.
## For precise Bayes factors, sampling at least 40,000 posterior
## samples is recommended.
## Warning in summarize_brms(pa_obj_log_pushing, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this
## was intended.
| exp(Est.) | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 106.16*** | 7.66 | [ 93.32, 122.54] | 1.000 | [0.94, 1.06] | 0.000 | >100 | Overwhelming Evidence | 1.003 | 1682 | 2179 |
| is_B | 146.96*** | 10.28 | [127.91, 167.84] | 1.000 | [0.94, 1.06] | 0.000 | >100 | Overwhelming Evidence | 1.001 | 2191 | 2869 |
| barriers_self_cw | 0.97 | 0.02 | [ 0.93, 1.00] | 0.970 | [0.94, 1.06] | 0.923 | 0.032 | Strong Evidence for Null | 1.001 | 9179 | 3053 |
| barriers_self_cb | 0.86 | 0.08 | [ 0.72, 1.05] | 0.928 | [0.94, 1.06] | 0.185 | 0.045 | Strong Evidence for Null | 1.001 | 2102 | 2746 |
| is_A:pushing_self_cw | 1.07 | 0.04 | [ 0.99, 1.15] | 0.958 | [0.94, 1.06] | 0.460 | 0.043 | Strong Evidence for Null | 1.001 | 4454 | 3162 |
| is_B:pushing_self_cw | 1.00 | 0.03 | [ 0.94, 1.07] | 0.508 | [0.94, 1.06] | 0.949 | 0.008 | Very Strong Evidence for Null | 1.000 | 4795 | 3193 |
| is_A:pushing_partner_cw | 0.98 | 0.02 | [ 0.94, 1.04] | 0.744 | [0.94, 1.06] | 0.963 | 0.007 | Very Strong Evidence for Null | 1.000 | 7078 | 2971 |
| is_B:pushing_partner_cw | 1.06 | 0.04 | [ 0.99, 1.14] | 0.949 | [0.94, 1.06] | 0.513 | 0.030 | Strong Evidence for Null | 1.000 | 5136 | 3496 |
| is_A:pushing_self_cb | 0.99 | 0.30 | [ 0.55, 1.87] | 0.511 | [0.94, 1.06] | 0.158 | 0.016 | Very Strong Evidence for Null | 1.002 | 1901 | 2345 |
| is_B:pushing_self_cb | 0.97 | 0.30 | [ 0.53, 1.85] | 0.538 | [0.94, 1.06] | 0.152 | 0.017 | Very Strong Evidence for Null | 1.002 | 1552 | 2479 |
| is_A:pushing_partner_cb | 1.19 | 0.38 | [ 0.62, 2.26] | 0.699 | [0.94, 1.06] | 0.136 | 0.019 | Very Strong Evidence for Null | 1.001 | 1502 | 2456 |
| is_B:pushing_partner_cb | 1.42 | 0.43 | [ 0.79, 2.60] | 0.875 | [0.94, 1.06] | 0.090 | 0.032 | Strong Evidence for Null | 1.002 | 2238 | 2707 |
| is_A:day | 1.01 | 0.06 | [ 0.88, 1.15] | 0.527 | [0.94, 1.06] | 0.659 | 0.005 | Very Strong Evidence for Null | 1.000 | 9192 | 3174 |
| is_B:day | 0.89 | 0.06 | [ 0.78, 1.01] | 0.960 | [0.94, 1.06] | 0.190 | 0.030 | Very Strong Evidence for Null | 1.004 | 9607 | 2369 |
| is_A:weartime_self_cw | 1.00** | 0.00 | [ 1.00, 1.00] | 0.997 | [0.94, 1.06] | 1.000 | 0.254 | Moderate Evidence for Null | 1.000 | 7285 | 3602 |
| is_B:weartime_self_cw | 1.00* | 0.00 | [ 1.00, 1.00] | 0.991 | [0.94, 1.06] | 1.000 | 0.102 | Moderate Evidence for Null | 1.000 | 8180 | 2831 |
| is_A:weartime_self_cb | 1.00 | 0.00 | [ 1.00, 1.00] | 0.677 | [0.94, 1.06] | 1.000 | 0.014 | Very Strong Evidence for Null | 1.003 | 1643 | 2402 |
| is_B:weartime_self_cb | 1.00 | 0.00 | [ 1.00, 1.00] | 0.875 | [0.94, 1.06] | 1.000 | 0.028 | Very Strong Evidence for Null | 1.000 | 2133 | 2900 |
| sd(is_A) | 0.34 | 0.05 | [0.26, 0.46] | NA | NA | NA | NA | NA | 1.004 | 1819 | 2307 |
| sd(is_B) | 0.34 | 0.05 | [0.26, 0.45] | NA | NA | NA | NA | NA | 1.001 | 1875 | 2972 |
| sd(is_A:pushing_self_cw) | 0.10 | 0.05 | [0.01, 0.22] | NA | NA | NA | NA | NA | 1.002 | 1096 | 1604 |
| sd(is_B:pushing_self_cw) | 0.09 | 0.04 | [0.02, 0.17] | NA | NA | NA | NA | NA | 1.001 | 1308 | 1555 |
| sd(is_A:pushing_partner_cw) | 0.03 | 0.03 | [0.00, 0.10] | NA | NA | NA | NA | NA | 1.000 | 2675 | 2695 |
| sd(is_B:pushing_partner_cw) | 0.10 | 0.04 | [0.01, 0.18] | NA | NA | NA | NA | NA | 1.001 | 1561 | 1168 |
| sigma | 0.52 | 0.01 | [0.50, 0.54] | NA | NA | NA | NA | NA | 1.001 | 6494 | 2778 |
## [1] 0 5
formula <- bf(
aff ~ 0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
is_A:persuasion_self_cb + is_B:persuasion_self_cb +
is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 5)", class = "b")
, brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = gaussian()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
mood_gauss_persuasion <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = gaussian(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("mood_gauss_persuasion", suffix))
)## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
check_brms(mood_gauss_persuasion, log_pp_check = FALSE)
DHARMa.check_brms.all(mood_gauss_persuasion, integer = FALSE)
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.15 [1.11, 1.19] 1.07 0.87
## is_B 1.04 [1.02, 1.09] 1.02 0.96
## is_A:persuasion_self_cw 1.06 [1.04, 1.11] 1.03 0.94
## is_B:persuasion_self_cw 1.10 [1.07, 1.14] 1.05 0.91
## is_A:persuasion_partner_cw 1.04 [1.02, 1.09] 1.02 0.96
## is_B:persuasion_partner_cw 2.66 [2.54, 2.80] 1.63 0.38
## is_A:persuasion_self_cb 2.69 [2.57, 2.83] 1.64 0.37
## Tolerance 95% CI
## [0.84, 0.90]
## [0.92, 0.98]
## [0.90, 0.96]
## [0.88, 0.94]
## [0.92, 0.98]
## [0.36, 0.39]
## [0.35, 0.39]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 47%
## normal 34%
## exponential 6%
##
## Predicted Distribution of Response
##
## Distribution Probability
## tweedie 41%
## beta-binomial 38%
## half-cauchy 6%
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: model.check
## outliers at both margin(s) = 38, observations = 3688, p-value =
## 1.198e-15
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.007301532 0.014115438
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.01030369
summary_mood <- summarize_brms(
mood_gauss_persuasion,
stats_to_report = stats_to_report,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = F) ## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
## Warning: Bayes factors might not be precise.
## For precise Bayes factors, sampling at least 40,000 posterior
## samples is recommended.
| Est. | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 3.61*** | 0.12 | [ 3.37, 3.84] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.002 | 453 | 899 |
| is_B | 3.87*** | 0.14 | [ 3.61, 4.16] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.010 | 446 | 1023 |
| is_A:persuasion_self_cw | 0.01 | 0.02 | [-0.04, 0.05] | 0.634 | [-0.11, 0.11] | 1.000 | 0.003 | Very Strong Evidence for Null | 1.000 | 5822 | 3143 |
| is_B:persuasion_self_cw | -0.01 | 0.02 | [-0.06, 0.03] | 0.708 | [-0.11, 0.11] | 1.000 | 0.003 | Very Strong Evidence for Null | 1.001 | 5420 | 2727 |
| is_A:persuasion_partner_cw | 0.05* | 0.02 | [ 0.00, 0.09] | 0.980 | [-0.11, 0.11] | 0.998 | 0.025 | Very Strong Evidence for Null | 1.000 | 5603 | 2712 |
| is_B:persuasion_partner_cw | 0.01 | 0.02 | [-0.04, 0.05] | 0.633 | [-0.11, 0.11] | 1.000 | 0.003 | Very Strong Evidence for Null | 1.000 | 5288 | 3058 |
| is_A:persuasion_self_cb | -0.27 | 0.32 | [-0.90, 0.39] | 0.797 | [-0.11, 0.11] | 0.190 | 0.023 | Very Strong Evidence for Null | 1.002 | 571 | 957 |
| is_B:persuasion_self_cb | 0.38 | 0.33 | [-0.31, 1.02] | 0.868 | [-0.11, 0.11] | 0.142 | 0.036 | Strong Evidence for Null | 1.003 | 672 | 1116 |
| is_A:persuasion_partner_cb | 0.39 | 0.27 | [-0.18, 0.92] | 0.917 | [-0.11, 0.11] | 0.120 | 0.040 | Strong Evidence for Null | 1.003 | 583 | 1141 |
| is_B:persuasion_partner_cb | -0.20 | 0.39 | [-0.97, 0.58] | 0.701 | [-0.11, 0.11] | 0.206 | 0.022 | Very Strong Evidence for Null | 1.005 | 672 | 1191 |
| is_A:day | 0.45*** | 0.07 | [ 0.31, 0.59] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 4830 | 2774 |
| is_B:day | 0.08 | 0.07 | [-0.05, 0.22] | 0.874 | [-0.11, 0.11] | 0.656 | 0.005 | Very Strong Evidence for Null | 1.000 | 5104 | 3328 |
| sd(is_A) | 0.70 | 0.09 | [0.55, 0.90] | NA | NA | NA | NA | NA | 1.011 | 605 | 1204 |
| sd(is_B) | 0.79 | 0.10 | [0.63, 1.03] | NA | NA | NA | NA | NA | 1.009 | 774 | 1225 |
| sd(is_A:pushing_self_cw) | 0.11 | 0.06 | [0.01, 0.24] | NA | NA | NA | NA | NA | 1.009 | 1039 | 1100 |
| sd(is_B:pushing_self_cw) | 0.06 | 0.05 | [0.00, 0.18] | NA | NA | NA | NA | NA | 1.002 | 1271 | 1356 |
| sd(is_A:pushing_partner_cw) | 0.13 | 0.05 | [0.03, 0.24] | NA | NA | NA | NA | NA | 1.007 | 732 | 301 |
| sd(is_B:pushing_partner_cw) | 0.12 | 0.06 | [0.01, 0.25] | NA | NA | NA | NA | NA | 1.003 | 1411 | 1092 |
| sigma | 0.87 | 0.01 | [0.85, 0.89] | NA | NA | NA | NA | NA | 1.000 | 5092 | 2948 |
formula <- bf(
aff ~ 0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw +
is_A:pressure_self_cb + is_B:pressure_self_cb +
is_A:pressure_partner_cb + is_B:pressure_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 5)", class = "b")
, brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = gaussian()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
mood_gauss_pressure <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = gaussian(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("mood_gauss_pressure", suffix))
)## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
check_brms(mood_gauss_pressure, log_pp_check = FALSE)
DHARMa.check_brms.all(mood_gauss_pressure, integer = FALSE)
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.20 [1.17, 1.25] 1.10 0.83
## is_B 1.01 [1.00, 1.12] 1.01 0.99
## is_A:pressure_self_cw 1.02 [1.00, 1.10] 1.01 0.98
## is_B:pressure_self_cw 1.03 [1.01, 1.09] 1.01 0.97
## is_A:pressure_partner_cw 1.00 [1.00, 3.28] 1.00 1.00
## Tolerance 95% CI
## [0.80, 0.86]
## [0.89, 1.00]
## [0.91, 1.00]
## [0.92, 0.99]
## [0.30, 1.00]
##
## Moderate Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_B:pressure_partner_cw 8.85 [8.37, 9.37] 2.98 0.11
## is_A:pressure_self_cb 8.07 [7.63, 8.53] 2.84 0.12
## Tolerance 95% CI
## [0.11, 0.12]
## [0.12, 0.13]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 47%
## normal 34%
## exponential 6%
##
## Predicted Distribution of Response
##
## Distribution Probability
## tweedie 41%
## beta-binomial 38%
## half-cauchy 6%
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: model.check
## outliers at both margin(s) = 38, observations = 3684, p-value =
## 1.158e-15
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.007309471 0.014130735
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.01031488
summary_mood <- summarize_brms(
mood_gauss_pressure,
stats_to_report = stats_to_report,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = F) ## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
## Warning: Bayes factors might not be precise.
## For precise Bayes factors, sampling at least 40,000 posterior
## samples is recommended.
| Est. | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 3.63*** | 0.14 | [ 3.35, 3.90] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.020 | 198 | 409 |
| is_B | 3.87*** | 0.14 | [ 3.58, 4.16] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.018 | 233 | 687 |
| is_A:pressure_self_cw | -0.06 | 0.06 | [-0.17, 0.06] | 0.844 | [-0.11, 0.11] | 0.841 | 0.005 | Very Strong Evidence for Null | 1.001 | 6955 | 3236 |
| is_B:pressure_self_cw | 0.01 | 0.06 | [-0.11, 0.11] | 0.542 | [-0.11, 0.11] | 0.960 | 0.003 | Very Strong Evidence for Null | 1.001 | 4185 | 3236 |
| is_A:pressure_partner_cw | 0.05 | 0.06 | [-0.06, 0.16] | 0.809 | [-0.11, 0.11] | 0.877 | 0.004 | Very Strong Evidence for Null | 1.002 | 3919 | 3108 |
| is_B:pressure_partner_cw | 0.04 | 0.06 | [-0.07, 0.15] | 0.752 | [-0.11, 0.11] | 0.900 | 0.004 | Very Strong Evidence for Null | 1.000 | 7136 | 2748 |
| is_A:pressure_self_cb | -0.18 | 0.67 | [-1.51, 1.16] | 0.602 | [-0.11, 0.11] | 0.125 | 0.017 | Very Strong Evidence for Null | 1.001 | 994 | 1789 |
| is_B:pressure_self_cb | 0.11 | 0.57 | [-1.02, 1.27] | 0.581 | [-0.11, 0.11] | 0.168 | 0.017 | Very Strong Evidence for Null | 1.007 | 1017 | 1699 |
| is_A:pressure_partner_cb | 0.23 | 0.52 | [-0.78, 1.28] | 0.675 | [-0.11, 0.11] | 0.157 | 0.015 | Very Strong Evidence for Null | 1.003 | 1025 | 1644 |
| is_B:pressure_partner_cb | -0.12 | 0.73 | [-1.59, 1.29] | 0.570 | [-0.11, 0.11] | 0.127 | 0.016 | Very Strong Evidence for Null | 1.005 | 977 | 1542 |
| is_A:day | 0.41*** | 0.07 | [ 0.27, 0.55] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 5190 | 2993 |
| is_B:day | 0.09 | 0.07 | [-0.05, 0.23] | 0.899 | [-0.11, 0.11] | 0.642 | 0.007 | Very Strong Evidence for Null | 1.001 | 5465 | 3359 |
| sd(is_A) | 0.71 | 0.09 | [0.56, 0.93] | NA | NA | NA | NA | NA | 1.009 | 432 | 926 |
| sd(is_B) | 0.81 | 0.10 | [0.63, 1.05] | NA | NA | NA | NA | NA | 1.005 | 610 | 888 |
| sd(is_A:pushing_self_cw) | 0.12 | 0.06 | [0.01, 0.25] | NA | NA | NA | NA | NA | 1.002 | 814 | 588 |
| sd(is_B:pushing_self_cw) | 0.06 | 0.05 | [0.00, 0.18] | NA | NA | NA | NA | NA | 1.003 | 1108 | 1004 |
| sd(is_A:pushing_partner_cw) | 0.13 | 0.05 | [0.04, 0.24] | NA | NA | NA | NA | NA | 1.002 | 1571 | 1131 |
| sd(is_B:pushing_partner_cw) | 0.12 | 0.06 | [0.02, 0.25] | NA | NA | NA | NA | NA | 1.002 | 1430 | 1166 |
| sigma | 0.87 | 0.01 | [0.85, 0.89] | NA | NA | NA | NA | NA | 1.001 | 5267 | 3080 |
formula <- bf(
aff ~ 0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw +
is_A:pushing_self_cb + is_B:pushing_self_cb +
is_A:pushing_partner_cb + is_B:pushing_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
barriers_self_cw + barriers_self_cb +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 5)", class = "b")
, brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = gaussian()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
mood_gauss_pushing <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = gaussian(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("mood_gauss_pushing", suffix))
)## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
check_brms(mood_gauss_pushing, log_pp_check = FALSE)
DHARMa.check_brms.all(mood_gauss_pushing, integer = FALSE)
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.38 [1.32, 1.46] 1.18 0.72
## is_B 1.01 [1.00, 1.39] 1.01 0.99
## barriers_self_cw 1.13 [1.09, 1.19] 1.06 0.89
## barriers_self_cb 1.09 [1.06, 1.15] 1.05 0.91
## is_A:pushing_self_cw 1.06 [1.03, 1.12] 1.03 0.95
## is_B:pushing_self_cw 1.07 [1.04, 1.14] 1.04 0.93
## is_A:pushing_partner_cw 1.12 [1.08, 1.18] 1.06 0.89
## is_B:pushing_partner_cw 1.51 [1.44, 1.60] 1.23 0.66
## is_A:pushing_self_cb 1.63 [1.55, 1.72] 1.28 0.61
## Tolerance 95% CI
## [0.68, 0.76]
## [0.72, 1.00]
## [0.84, 0.92]
## [0.87, 0.95]
## [0.89, 0.97]
## [0.88, 0.96]
## [0.85, 0.93]
## [0.63, 0.70]
## [0.58, 0.65]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## cauchy 41%
## normal 41%
## exponential 6%
##
## Predicted Distribution of Response
##
## Distribution Probability
## tweedie 41%
## beta-binomial 38%
## half-cauchy 6%
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: model.check
## outliers at both margin(s) = 15, observations = 1776, p-value =
## 4.841e-06
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.004734615 0.013892099
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.008445946
summary_mood <- summarize_brms(
mood_gauss_pushing,
stats_to_report = stats_to_report,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = F) ## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
## Warning: Bayes factors might not be precise.
## For precise Bayes factors, sampling at least 40,000 posterior
## samples is recommended.
| Est. | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 3.65*** | 0.12 | [ 3.38, 3.89] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.002 | 826 | 1559 |
| is_B | 3.84*** | 0.15 | [ 3.54, 4.13] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.006 | 918 | 1668 |
| barriers_self_cw | -0.21*** | 0.03 | [-0.26, -0.16] | 1.000 | [-0.11, 0.11] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 8770 | 3006 |
| barriers_self_cb | -0.35 | 0.22 | [-0.80, 0.12] | 0.936 | [-0.11, 0.11] | 0.122 | 0.050 | Strong Evidence for Null | 1.005 | 789 | 1517 |
| is_A:pushing_self_cw | 0.09 | 0.05 | [-0.01, 0.18] | 0.964 | [-0.11, 0.11] | 0.671 | 0.030 | Strong Evidence for Null | 1.002 | 4384 | 3110 |
| is_B:pushing_self_cw | -0.01 | 0.04 | [-0.09, 0.06] | 0.644 | [-0.11, 0.11] | 0.990 | 0.005 | Very Strong Evidence for Null | 1.000 | 5420 | 3240 |
| is_A:pushing_partner_cw | 0.09 | 0.05 | [ 0.00, 0.18] | 0.973 | [-0.11, 0.11] | 0.674 | 0.040 | Strong Evidence for Null | 1.000 | 3733 | 2884 |
| is_B:pushing_partner_cw | 0.12* | 0.05 | [ 0.01, 0.23] | 0.984 | [-0.11, 0.11] | 0.387 | 0.077 | Strong Evidence for Null | 1.000 | 3298 | 2882 |
| is_A:pushing_self_cb | -0.17 | 0.49 | [-1.19, 0.79] | 0.637 | [-0.11, 0.11] | 0.165 | 0.014 | Very Strong Evidence for Null | 1.003 | 1189 | 1590 |
| is_B:pushing_self_cb | 0.82 | 0.68 | [-0.53, 2.25] | 0.891 | [-0.11, 0.11] | 0.065 | 0.043 | Strong Evidence for Null | 1.001 | 996 | 1543 |
| is_A:pushing_partner_cb | 0.78 | 0.56 | [-0.35, 1.92] | 0.916 | [-0.11, 0.11] | 0.054 | 0.040 | Strong Evidence for Null | 1.004 | 1020 | 1462 |
| is_B:pushing_partner_cb | -0.33 | 0.63 | [-1.58, 1.01] | 0.702 | [-0.11, 0.11] | 0.110 | 0.023 | Very Strong Evidence for Null | 1.002 | 1156 | 1592 |
| is_A:day | 0.35*** | 0.10 | [ 0.16, 0.54] | 1.000 | [-0.11, 0.11] | 0.007 | 1.512 | Weak Evidence | 1.001 | 7719 | 2951 |
| is_B:day | 0.08 | 0.10 | [-0.11, 0.28] | 0.806 | [-0.11, 0.11] | 0.589 | 0.006 | Very Strong Evidence for Null | 1.004 | 6465 | 2807 |
| sd(is_A) | 0.66 | 0.09 | [0.52, 0.87] | NA | NA | NA | NA | NA | 1.005 | 959 | 1597 |
| sd(is_B) | 0.80 | 0.10 | [0.63, 1.05] | NA | NA | NA | NA | NA | 1.002 | 1041 | 1715 |
| sd(is_A:pushing_self_cw) | 0.10 | 0.07 | [0.01, 0.24] | NA | NA | NA | NA | NA | 1.003 | 1278 | 1291 |
| sd(is_B:pushing_self_cw) | 0.07 | 0.06 | [0.00, 0.20] | NA | NA | NA | NA | NA | 1.000 | 1483 | 1803 |
| sd(is_A:pushing_partner_cw) | 0.12 | 0.05 | [0.02, 0.26] | NA | NA | NA | NA | NA | 1.001 | 1844 | 1640 |
| sd(is_B:pushing_partner_cw) | 0.13 | 0.07 | [0.01, 0.30] | NA | NA | NA | NA | NA | 1.002 | 1520 | 1406 |
| sigma | 0.84 | 0.01 | [0.81, 0.86] | NA | NA | NA | NA | NA | 1.001 | 5401 | 3024 |
## [1] 0 5
introduce_binary_reactance <- function(data) {
data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
return(data)
}
df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
for (i in seq_along(implist)) {
implist[[i]] <- introduce_binary_reactance(implist[[i]])
}
}formula <- bf(
is_reactance ~ 0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
is_A:persuasion_self_cb + is_B:persuasion_self_cb +
is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:persuasion_self_cw + is_B:persuasion_self_cw +
is_A:persuasion_partner_cw + is_B:persuasion_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 10)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 10)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = bernoulli()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
is_reactance_persuasion <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::bernoulli(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("is_reactance_persuasion", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
check_brms(is_reactance_persuasion)
DHARMa.check_brms.all(is_reactance_persuasion, integer = FALSE)
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.80 [1.69, 1.93] 1.34 0.55
## is_B 1.37 [1.30, 1.46] 1.17 0.73
## is_A:persuasion_self_cw 1.17 [1.11, 1.24] 1.08 0.86
## is_B:persuasion_self_cw 1.05 [1.02, 1.14] 1.02 0.95
## is_A:persuasion_partner_cw 1.06 [1.02, 1.14] 1.03 0.94
## is_B:persuasion_partner_cw 1.61 [1.51, 1.72] 1.27 0.62
## is_A:persuasion_self_cb 2.00 [1.87, 2.15] 1.41 0.50
## Tolerance 95% CI
## [0.52, 0.59]
## [0.68, 0.77]
## [0.81, 0.90]
## [0.88, 0.98]
## [0.88, 0.98]
## [0.58, 0.66]
## [0.47, 0.54]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## normal 44%
## cauchy 12%
## gamma 12%
##
## Predicted Distribution of Response
##
## Distribution Probability
## bernoulli 81%
## beta-binomial 12%
## binomial 6%
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 9 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: model.check
## outliers at both margin(s) = 2, observations = 756, p-value =
## 0.6663
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.0003205434 0.0095234999
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.002645503
summary_is_reactance <- summarize_brms(
is_reactance_persuasion,
stats_to_report = stats_to_report,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
## Warning: Bayes factors might not be precise.
## For precise Bayes factors, sampling at least 40,000 posterior
## samples is recommended.
| OR | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 0.18*** | 0.09 | [0.06, 0.42] | 1.000 | [0.83, 1.20] | 0.000 | 14.783 | Strong Evidence | 1.000 | 2614 | 2734 |
| is_B | 0.29** | 0.12 | [0.12, 0.64] | 0.999 | [0.83, 1.20] | 0.005 | 1.740 | Weak Evidence | 1.001 | 3369 | 3196 |
| is_A:persuasion_self_cw | 0.90 | 0.13 | [0.64, 1.17] | 0.769 | [0.83, 1.20] | 0.681 | 0.066 | Strong Evidence for Null | 1.001 | 2236 | 2875 |
| is_B:persuasion_self_cw | 0.78 | 0.11 | [0.54, 1.02] | 0.965 | [0.83, 1.20] | 0.308 | 0.309 | Weak Evidence for Null | 1.000 | 3027 | 2084 |
| is_A:persuasion_partner_cw | 0.88 | 0.17 | [0.59, 1.28] | 0.740 | [0.83, 1.20] | 0.573 | 0.084 | Strong Evidence for Null | 1.001 | 3281 | 2419 |
| is_B:persuasion_partner_cw | 1.26 | 0.24 | [0.85, 1.93] | 0.889 | [0.83, 1.20] | 0.370 | 0.154 | Moderate Evidence for Null | 1.001 | 2604 | 2514 |
| is_A:persuasion_self_cb | 5.18 | 4.43 | [0.95, 36.77] | 0.970 | [0.83, 1.20] | 0.029 | 0.520 | Weak Evidence for Null | 1.001 | 1970 | 2558 |
| is_B:persuasion_self_cb | 2.81 | 2.20 | [0.58, 14.93] | 0.905 | [0.83, 1.20] | 0.074 | 0.237 | Moderate Evidence for Null | 1.002 | 1937 | 2597 |
| is_A:persuasion_partner_cb | 1.54 | 1.10 | [0.39, 8.13] | 0.734 | [0.83, 1.20] | 0.175 | 0.106 | Moderate Evidence for Null | 1.002 | 1658 | 1952 |
| is_B:persuasion_partner_cb | 2.01 | 1.82 | [0.33, 14.67] | 0.784 | [0.83, 1.20] | 0.120 | 0.113 | Moderate Evidence for Null | 1.000 | 2156 | 2685 |
| is_A:day | 2.54 | 1.40 | [0.88, 7.62] | 0.959 | [0.83, 1.20] | 0.061 | 0.178 | Moderate Evidence for Null | 1.001 | 5826 | 3180 |
| is_B:day | 1.04 | 0.58 | [0.34, 3.16] | 0.528 | [0.83, 1.20] | 0.248 | 0.050 | Strong Evidence for Null | 1.003 | 5441 | 2933 |
| sd(is_A) | 1.37 | 0.42 | [0.67, 2.43] | NA | NA | NA | NA | NA | 1.006 | 1215 | 1900 |
| sd(is_B) | 1.37 | 0.33 | [0.84, 2.17] | NA | NA | NA | NA | NA | 1.000 | 2191 | 2654 |
| sd(is_A:persuasion_self_cw) | 0.30 | 0.18 | [0.02, 0.71] | NA | NA | NA | NA | NA | 1.003 | 794 | 963 |
| sd(is_B:persuasion_self_cw) | 0.35 | 0.26 | [0.02, 0.95] | NA | NA | NA | NA | NA | 1.002 | 797 | 1734 |
| sd(is_A:persuasion_partner_cw) | 0.48 | 0.27 | [0.05, 1.16] | NA | NA | NA | NA | NA | 1.002 | 1413 | 1509 |
| sd(is_B:persuasion_partner_cw) | 0.64 | 0.28 | [0.17, 1.37] | NA | NA | NA | NA | NA | 1.003 | 1425 | 1048 |
formula <- bf(
is_reactance ~ 0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw +
is_A:pressure_self_cb + is_B:pressure_self_cb +
is_A:pressure_partner_cb + is_B:pressure_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pressure_self_cw + is_B:pressure_self_cw +
is_A:pressure_partner_cw + is_B:pressure_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 10)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 10)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = bernoulli()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
is_reactance_pressure <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::bernoulli(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("is_reactance_pressure", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
check_brms(is_reactance_pressure)
DHARMa.check_brms.all(is_reactance_pressure, integer = FALSE)
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.82 [1.70, 1.95] 1.35 0.55
## is_B 1.03 [1.01, 1.15] 1.02 0.97
## is_A:pressure_self_cw 1.02 [1.00, 1.29] 1.01 0.98
## is_B:pressure_self_cw 1.01 [1.00, 1.37] 1.01 0.99
## is_A:pressure_partner_cw 1.01 [1.00, 2.53] 1.00 0.99
## is_B:pressure_partner_cw 2.48 [2.30, 2.67] 1.57 0.40
## is_A:pressure_self_cb 2.93 [2.71, 3.17] 1.71 0.34
## Tolerance 95% CI
## [0.51, 0.59]
## [0.87, 0.99]
## [0.78, 1.00]
## [0.73, 1.00]
## [0.39, 1.00]
## [0.37, 0.43]
## [0.32, 0.37]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## normal 34%
## beta 16%
## cauchy 12%
##
## Predicted Distribution of Response
##
## Distribution Probability
## bernoulli 81%
## beta-binomial 12%
## binomial 6%
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 26 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: model.check
## outliers at both margin(s) = 1, observations = 756, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.0000334886 0.0073476538
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.001322751
summary_is_reactance <- summarize_brms(
is_reactance_pressure,
stats_to_report = stats_to_report,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
## Warning: Bayes factors might not be precise.
## For precise Bayes factors, sampling at least 40,000 posterior
## samples is recommended.
| OR | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 0.15*** | 0.06 | [0.06, 0.32] | 1.000 | [0.83, 1.20] | 0.000 | 87.472 | Very Strong Evidence | 1.001 | 2317 | 2099 |
| is_B | 0.21*** | 0.07 | [0.10, 0.39] | 1.000 | [0.83, 1.20] | 0.000 | 95.261 | Very Strong Evidence | 1.000 | 3422 | 2844 |
| is_A:pressure_self_cw | 2.39 | 1.22 | [0.79, 10.38] | 0.951 | [0.83, 1.20] | 0.053 | 0.690 | Weak Evidence for Null | 1.001 | 1952 | 1596 |
| is_B:pressure_self_cw | 1.87 | 0.65 | [0.82, 4.46] | 0.943 | [0.83, 1.20] | 0.087 | 0.354 | Weak Evidence for Null | 1.000 | 2479 | 1818 |
| is_A:pressure_partner_cw | 0.87 | 0.66 | [0.09, 3.64] | 0.576 | [0.83, 1.20] | 0.202 | 0.127 | Moderate Evidence for Null | 1.002 | 2272 | 2073 |
| is_B:pressure_partner_cw | 1.55 | 0.76 | [0.42, 5.01] | 0.807 | [0.83, 1.20] | 0.176 | 0.129 | Moderate Evidence for Null | 1.001 | 3059 | 2828 |
| is_A:pressure_self_cb | 17.11 | 29.38 | [0.80, 1025.83] | 0.967 | [0.83, 1.20] | 0.015 | 0.601 | Weak Evidence for Null | 1.002 | 2273 | 1848 |
| is_B:pressure_self_cb | 22.64** | 32.84 | [1.90, 815.16] | 0.996 | [0.83, 1.20] | 0.005 | 2.398 | Weak Evidence | 1.000 | 2655 | 2128 |
| is_A:pressure_partner_cb | 1.05 | 1.45 | [0.05, 16.96] | 0.514 | [0.83, 1.20] | 0.106 | 0.129 | Moderate Evidence for Null | 1.002 | 2319 | 1259 |
| is_B:pressure_partner_cb | 0.36 | 0.64 | [0.00, 8.17] | 0.725 | [0.83, 1.20] | 0.073 | 0.104 | Moderate Evidence for Null | 1.001 | 3171 | 2356 |
| is_A:day | 2.26 | 1.17 | [0.81, 6.51] | 0.943 | [0.83, 1.20] | 0.074 | 0.147 | Moderate Evidence for Null | 1.000 | 5731 | 2868 |
| is_B:day | 1.33 | 0.70 | [0.47, 3.72] | 0.705 | [0.83, 1.20] | 0.226 | 0.052 | Strong Evidence for Null | 1.002 | 4981 | 2411 |
| sd(is_A) | 1.46 | 0.39 | [0.83, 2.40] | NA | NA | NA | NA | NA | 1.001 | 1456 | 2124 |
| sd(is_B) | 1.14 | 0.28 | [0.66, 1.84] | NA | NA | NA | NA | NA | 1.000 | 1388 | 1435 |
| sd(is_A:pressure_self_cw) | 1.28 | 1.08 | [0.05, 3.99] | NA | NA | NA | NA | NA | 1.003 | 872 | 1645 |
| sd(is_B:pressure_self_cw) | 1.02 | 0.57 | [0.17, 2.51] | NA | NA | NA | NA | NA | 1.002 | 1315 | 1027 |
| sd(is_A:pressure_partner_cw) | 1.81 | 1.16 | [0.15, 4.45] | NA | NA | NA | NA | NA | 1.001 | 1296 | 1321 |
| sd(is_B:pressure_partner_cw) | 0.71 | 0.64 | [0.03, 2.84] | NA | NA | NA | NA | NA | 1.001 | 2176 | 1877 |
formula <- bf(
is_reactance ~ 0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw +
is_A:pushing_self_cb + is_B:pushing_self_cb +
is_A:pushing_partner_cb + is_B:pushing_partner_cb +
#is_A:plan_self + is_B:plan_self +
#is_A:plan_partner + is_B:plan_partner +
barriers_self_cw + barriers_self_cb +
is_A:day + is_B:day +
# Random effects
(0 +
is_A + is_B + # Intercepts
is_A:pushing_self_cw + is_B:pushing_self_cw +
is_A:pushing_partner_cw + is_B:pushing_partner_cw |dd| coupleID)
, decomp = 'QR'
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 10)", class = "b", coef = 'is_A')
, brms::set_prior("normal(0, 10)", class = "b", coef = 'is_B')
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = bernoulli()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
is_reactance_pushing <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::bernoulli(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("is_reactance_pushing", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.
if (check_models) {
check_brms(is_reactance_pushing)
DHARMa.check_brms.all(is_reactance_pushing, integer = FALSE)
}## Model has no intercept. VIFs may not be sensible.
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## is_A 1.89 [1.78, 2.01] 1.37 0.53
## is_B 1.07 [1.03, 1.13] 1.03 0.94
## barriers_self_cw 1.11 [1.07, 1.17] 1.05 0.90
## barriers_self_cb 1.24 [1.18, 1.31] 1.11 0.81
## is_A:pushing_self_cw 1.12 [1.08, 1.18] 1.06 0.89
## is_B:pushing_self_cw 1.02 [1.00, 1.21] 1.01 0.98
## is_A:pushing_partner_cw 1.05 [1.02, 1.12] 1.03 0.95
## is_B:pushing_partner_cw 1.47 [1.40, 1.56] 1.21 0.68
## is_A:pushing_self_cb 1.37 [1.31, 1.45] 1.17 0.73
## Tolerance 95% CI
## [0.50, 0.56]
## [0.88, 0.97]
## [0.86, 0.94]
## [0.77, 0.84]
## [0.85, 0.93]
## [0.82, 1.00]
## [0.89, 0.98]
## [0.64, 0.71]
## [0.69, 0.77]
## # Distribution of Model Family
##
## Predicted Distribution of Residuals
##
## Distribution Probability
## normal 44%
## cauchy 12%
## gamma 9%
##
## Predicted Distribution of Response
##
## Distribution Probability
## bernoulli 75%
## beta-binomial 16%
## binomial 9%
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 8 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: model.check
## outliers at both margin(s) = 2, observations = 486, p-value =
## 0.2536
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.0004987621 0.0147859211
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.004115226
summary_is_reactance <- summarize_brms(
is_reactance_pushing,
stats_to_report = stats_to_report,
#model_rows_fixed = model_rows_fixed,
#model_rows_random = model_rows_random,
#model_rownames_fixed = model_rownames_fixed,
#model_rownames_random = model_rownames_random,
exponentiate = T) ## Model has no intercept. VIFs may not be sensible.
## Sampling priors, please wait...
## Warning: Bayes factors might not be precise.
## For precise Bayes factors, sampling at least 40,000 posterior
## samples is recommended.
| OR | SE | 95% CI | pd | ROPE | inside ROPE | BF | BF_Evidence | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| is_A | 0.13*** | 0.07 | [0.04, 0.33] | 1.000 | [0.83, 1.20] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 2447 | 2234 |
| is_B | 0.16*** | 0.07 | [0.06, 0.37] | 1.000 | [0.83, 1.20] | 0.000 | >100 | Overwhelming Evidence | 1.000 | 3196 | 2361 |
| barriers_self_cw | 0.78 | 0.14 | [0.54, 1.10] | 0.919 | [0.83, 1.20] | 0.342 | 0.139 | Moderate Evidence for Null | 1.000 | 5070 | 2919 |
| barriers_self_cb | 2.48 | 1.81 | [0.62, 10.60] | 0.900 | [0.83, 1.20] | 0.096 | 0.211 | Moderate Evidence for Null | 1.001 | 2695 | 2661 |
| is_A:pushing_self_cw | 1.35 | 0.25 | [0.93, 2.02] | 0.944 | [0.83, 1.20] | 0.249 | 0.232 | Moderate Evidence for Null | 1.000 | 4243 | 2725 |
| is_B:pushing_self_cw | 1.44 | 0.30 | [0.97, 2.27] | 0.964 | [0.83, 1.20] | 0.157 | 0.463 | Weak Evidence for Null | 1.001 | 3100 | 2454 |
| is_A:pushing_partner_cw | 0.77 | 0.24 | [0.31, 1.33] | 0.814 | [0.83, 1.20] | 0.342 | 0.119 | Moderate Evidence for Null | 1.003 | 2255 | 1751 |
| is_B:pushing_partner_cw | 1.16 | 0.24 | [0.76, 1.74] | 0.761 | [0.83, 1.20] | 0.500 | 0.080 | Strong Evidence for Null | 1.001 | 4260 | 3039 |
| is_A:pushing_self_cb | 29.91* | 47.79 | [1.80, 1723.06] | 0.992 | [0.83, 1.20] | 0.008 | 1.715 | Weak Evidence | 1.002 | 1684 | 1694 |
| is_B:pushing_self_cb | 2.52 | 4.01 | [0.11, 74.90] | 0.725 | [0.83, 1.20] | 0.080 | 0.125 | Moderate Evidence for Null | 1.001 | 2147 | 2542 |
| is_A:pushing_partner_cb | 0.08 | 0.16 | [0.00, 3.84] | 0.898 | [0.83, 1.20] | 0.029 | 0.237 | Moderate Evidence for Null | 1.001 | 2349 | 2793 |
| is_B:pushing_partner_cb | 1.20 | 1.87 | [0.04, 27.96] | 0.545 | [0.83, 1.20] | 0.092 | 0.107 | Moderate Evidence for Null | 1.001 | 2505 | 2772 |
| is_A:day | 1.31 | 0.93 | [0.32, 5.22] | 0.658 | [0.83, 1.20] | 0.192 | 0.055 | Strong Evidence for Null | 1.002 | 4646 | 2904 |
| is_B:day | 1.38 | 0.92 | [0.37, 5.25] | 0.685 | [0.83, 1.20] | 0.199 | 0.065 | Strong Evidence for Null | 1.001 | 6392 | 2678 |
| sd(is_A) | 1.45 | 0.49 | [0.61, 2.66] | NA | NA | NA | NA | NA | 1.000 | 1365 | 2100 |
| sd(is_B) | 1.30 | 0.40 | [0.65, 2.26] | NA | NA | NA | NA | NA | 1.001 | 1836 | 2125 |
| sd(is_A:pushing_self_cw) | 0.25 | 0.22 | [0.01, 0.86] | NA | NA | NA | NA | NA | 1.003 | 1342 | 1956 |
| sd(is_B:pushing_self_cw) | 0.57 | 0.31 | [0.07, 1.39] | NA | NA | NA | NA | NA | 1.008 | 1046 | 1246 |
| sd(is_A:pushing_partner_cw) | 0.64 | 0.49 | [0.04, 1.92] | NA | NA | NA | NA | NA | 1.005 | 1150 | 1424 |
| sd(is_B:pushing_partner_cw) | 0.29 | 0.25 | [0.01, 1.06] | NA | NA | NA | NA | NA | 1.001 | 1769 | 2283 |
Analyses were conducted using the R Statistical language (version 4.4.2; R Core Team, 2024) on Windows 11 x64 (build 26100)